# Estimate models using alternative time windows (±5, 10, 20, 25, 30 days).
# See Appendix K.




######################################
# loop through different time ranges #
######################################

# function to fit stan models with different time windows
estTime <- function(dayrange){
  library(brms)
  library(ggplot2)
  
  # load data
  load("PSRM Replication Files/PresidentialElectionsPosts.RData")
  
  # subset data by date range
  day_range <- dayrange
  
  filename <- paste0("Stan_", dayrange, "Days.RData")
  
  sub <- president_df[which(president_df$daysinceelection >= -1*day_range
                            & president_df$daysinceelection <= day_range),]
  
  winners <- sub[which(sub$win == 1),]
  losers <- sub[which(sub$win != 1),]
  
  # print data summary
  print(day_range)
  print(filename)
  print(c(nrow(winners), nrow(losers)))
  
  # winners + love
  winner_love2 <- brm(loveprop ~ post*populist_dummy + post*polarization
                      + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                      + Xl + Xr 
                      + (1 + Xl + Xr | party_election),
                      data=winners,
                      warmup=1000, iter=2000, chains=3, seed=1234,
                      silent = 2, refresh = 0)
  
  # winners + angry
  winner_angry2 <- brm(angryprop ~ post*populist_dummy + post*polarization
                       + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                       + Xl + Xr 
                       + (1 + Xl + Xr | party_election),
                       data=winners,
                       warmup=1000, iter=2000, chains=3, seed=1234,
                       silent = 2, refresh = 0)
  
  # losers + love
  loser_love2 <- brm(loveprop ~ post*populist_dummy + post*polarization
                     + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                     + Xl + Xr 
                     + (1 + Xl + Xr | party_election),
                     data=losers,
                     warmup=1000, iter=2000, chains=3, seed=1234,
                     silent = 2, refresh = 0)

  # losers + angry
  loser_angry2 <- brm(angryprop ~ post*populist_dummy + post*polarization
                      + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                      + Xl + Xr 
                      + (1 + Xl + Xr | party_election),
                      data=losers,
                      warmup=1000, iter=2000, chains=3, seed=1234,
                      silent = 2, refresh = 0)

  save(winner_love2,
       winner_angry2,
       loser_love2,
       loser_angry2,
       file=filename)
}

estTime(5) # 1531 1108
estTime(10) # 2569 1969
estTime(20) # 4664 3658
estTime(25) # 5618 4526
estTime(30) # 6577 5439




###########
# summary #
###########

# check convergence
load("PSRM Replication Files/Stan_5Days.RData")
length(ls())
max(summary(winner_love2)[[14]][,5])
max(summary(winner_angry2)[[14]][,5])
max(summary(loser_love2)[[14]][,5])
max(summary(loser_angry2)[[14]][,5])
rm(list=ls())

load("PSRM Replication Files/Stan_10Days.RData")
length(ls())
max(summary(winner_love2)[[14]][,5])
max(summary(winner_angry2)[[14]][,5])
max(summary(loser_love2)[[14]][,5])
max(summary(loser_angry2)[[14]][,5])
rm(list=ls())

load("PSRM Replication Files/Stan_20Days.RData")
length(ls())
max(summary(winner_love2)[[14]][,5])
max(summary(winner_angry2)[[14]][,5])
max(summary(loser_love2)[[14]][,5])
max(summary(loser_angry2)[[14]][,5])
rm(list=ls())

load("PSRM Replication Files/Stan_25Days.RData")
length(ls())
max(summary(winner_love2)[[14]][,5])
max(summary(winner_angry2)[[14]][,5])
max(summary(loser_love2)[[14]][,5])
max(summary(loser_angry2)[[14]][,5])
rm(list=ls())

load("PSRM Replication Files/Stan_30Days.RData")
length(ls())
max(summary(winner_love2)[[14]][,5])
max(summary(winner_angry2)[[14]][,5])
max(summary(loser_love2)[[14]][,5])
max(summary(loser_angry2)[[14]][,5])
rm(list=ls())


# load stan model outputs with different time windows
folder <- "PSRM Replication Files"
(files <- list.files(folder, full.names = TRUE))
(files <- files[grep("Days", files)])

citab <- data.frame(NULL)

timerange <- c(15,10,20,25,30,5)

# get credible intervals
getCI <- function(model, intr){
  coeftab <- summary(model)[[14]][,c(1,3,4)]
  coeftab <- coeftab[intr,]
  return(coeftab)
}

for(i in 1:6){
  load(files[i])
  
  dtemp <- data.frame(rbind(getCI(winner_love2, 12), getCI(winner_angry2, 12), 
                            getCI(winner_love2, 13), getCI(winner_angry2, 13),
                            getCI(loser_love2, 12), getCI(loser_angry2, 12),
                            getCI(loser_love2, 13), getCI(loser_angry2, 13)))
  
  colnames(dtemp) <- c("est", "lwr", "upr")
  
  dtemp$type <- c("Winner + Love + Populist Involvement", "Winner + Angry + Populist Involvement", 
                  "Winner + Love + Polarization", "Winner + Angry + Polarization",
                  "Loser + Love + Populist Involvement", "Loser + Angry + Populist Involvement",
                  "Loser + Love + Polarization", "Loser + Angry + Polarization")
  
  dtemp$timerange <- timerange[i]
  
  citab <- rbind(citab, dtemp)
}

table(citab$timerange)

citab$color <- "black"
citab$color[citab$timerange == 15] <- "red"

est_types <- unique(citab$type)
(est_types <- est_types[c(1,2,5,6,3,4,7,8)])

# figure K.1
#pdf("timewindow.pdf", width=10, height=14)
par(mfrow=c(4,2))
for(i in est_types){
  sub <- citab[which(citab$type == i),]
  sub <- sub[order(sub$timerange),]
  
  plot(nrow(sub):1, sub$est, type="l", col="gray70",
       main=i,
       ylab="Posterior Estimate of the Interaction Term",
       xaxt="n", xlab="Day Range (±)",
       xlim=c(0.8,6.2),
       ylim=c(min(sub$lwr)-0.2, max(sub$upr)+0.2))
  abline(h=0, lty=2)
  axis(1, at=1:6, label=c(30, 25, 20, 15, 10, 5))
  points(nrow(sub):1, sub$est, pch=20, cex=2, col=sub$color)
  segments(nrow(sub):1, sub$lwr, nrow(sub):1, sub$upr, lwd=1.5, col=sub$color)
}
#dev.off()
